# -*- coding: utf-8 -*-
# Code sous licence CeCill
# http://www.cecill.info/licences.en.html
# This is a license compatible with the GNU GPL license 
# and that conforms to the principles of French law.
# Authors: 
#       Nicolas Bruyant, CNRS,LNCMI
"""
TDO experiment signal treatment
=============

Functions candidates for frequency trancking and signal characterization

Functions index:
    ZeroCrossingTDO: Recover instantaneous TDO signal by detecting zero crossing
            
    

TODO:
----
    

"""
__version__ = "Version 0.1, SVN-$Revision: 18 $"
"""
Revision allow a good tracking of data processing  errors
Refers to SVN revisions: DO NOT edit manually
"""
import numpy as np

#from tptsig.wavfile import *
#from tptsig.nmf import *
from scipy import *
from scipy.signal import *
from scipy.optimize import *
import pylab as plt
from importfilters import *
import matplotlib.mlab as mlab
import time






def ZeroCrossingTDO(y):
    """
    Recover instantaneous TDO signal by detecting zero crossing
    
    Parameters
    ==========
        Filename of a Gagescope file,
        TODO:
        Sample rate of the datas,
        Number of column of the datafile
        
    Returns
    =======
        Plots of the recovered signal
        Numpy array containing datas
        
    Examples
    ========  
    >>> filename="/Users/Nicolas/Documents/Projets/TDO/post traitement/tracking freq/s120725d.sig"
    >>> temp=loadfile(inputfile,10E6,4)
        
    """
    #TODO 
    
    # Detect zero crossing indexes
    # TOTEST: old zc = np.where(np.diff(np.sign(y))=1)[0]
    zc = np.where(np.diff(np.sign(y))>0)[0]
    
    # raw intervals datas
    #ints=zc[range(1,len(zc),1)]-c[range(0,len(zc)-1,1)]
    
    # calculate the zero crossing point
    by=zc+(abs(y[zc]))/(abs(y[zc]-y[zc+1])+0.0)
    
    
    # Calculate intervalls
    bya=by[:len(by)-1]
    byb=by[1:]
    intsb=byb-bya
    
    #average positives and negatives intervals (useless if an entire period is taken)
    #intsb2=(intsb[range(0,len(intsb)-1,2)]+intsb[range(1,len(intsb),2)])
    figure()
    plot(intsb)
    
    
    # Data filtering
    cutoff = 0.1
    filterorder = 5
    
    b,a = filter_design.butter(filterorder,cutoff,btype='low') 
    fints=filtfilt(b,a,intsb)
    plot(fints)
    
    return fints,intsb


def FFTZeropadTDO(sig):
    # To reduce signals size
    # y=y[:5002]
    (Pxx, freqs, bins)=mlab.specgram(sig, NFFT=2048,Fs=2,
               detrend=mlab.detrend_none, window=mlab.window_hanning,noverlap=1024,
               pad_to=2048*32,sides='default',scale_by_freq=None)
    
    #Freq=map(lambda x: freqs[np.argmax(x)],Pxx)
    #Amp=map(lambda x: np.max(x),Pxx)
    Amp=np.max(Pxx,axis=0)
    Freq=freqs[np.argmax(Pxx,axis=0)]
    return Amp,Freq

    #    # manual spectrogramm
    #    Nbslices=10000
    #    # slices
    #    SlicesDim=(Nbslices,size(y0)/Nbslices)
    #    y0.resize(SlicesDim)
    #    
    #    zeropad=2**12
    #    
    #    s1=size(np.fft.rfft(y0[1,:],n=zeropad))
    #        
    #    stfft=sp.zeros((s1,Nbslices),dtype=complex,)
    #        
    #    for i in range(Nbslices):
    #        stfft[:,i]=np.fft.rfft(y0[i,:],n=zeropad)
    #        
    #    power=log(abs(stfft))**2    
    #        
    #    plt.imshow(power[:300,:5000],origin='lower',aspect='auto')
    
def FFTGaussianfit(sig):
    (Pxx, freqs, bins)=mlab.specgram(sig, NFFT=1024,Fs=2,
               detrend=mlab.detrend_none, window=mlab.window_hanning,noverlap=256,
               pad_to=2048*2,sides='default',scale_by_freq=None)
    #(Pxx, freqs, bins, im)=plt.specgram(sig, NFFT=4096, noverlap=1024,pad_to=4096*1)            
    #figure(num=1)
    #plot(Pxx[:10,:1000].T)
    #figure(num=2)
    #imshow(log(abs(Pxx[:600,:])),origin='lower',aspect='auto')
    Amp=np.empty((size(Pxx,1),))
    Freq=np.empty((size(Pxx,1),))
    FreqMax=np.empty((size(Pxx,1),))
    Sigma=np.empty((size(Pxx,1),))
    for i in range(0,size(Pxx,1),1):
        #print i
        y=Pxx[:,i]
        #plot(freqs,log(y))
        popt, pcov,Max = MaxGaussianFit(y,freqs, debug = 0)
        # append to data list
        Amp[i]=popt[0]
        Freq[i]=popt[1]
        FreqMax[i]=Max
        Sigma[i]=popt[2]
        
    return Amp,Freq,FreqMax,Sigma    
        


def MaxGaussianFit(y,x, debug = 0):        
        indmax=np.argmax(y)
        ymax=np.max(y)
        
        # largeur du fit
        Indexs=arange(indmax-5,indmax+5)
        
        func = lambda x, k, tau, m: k* exp(-(x - tau) ** 2/m**2) 
        tau = x[indmax]
        # get a first approximation of peak amplitude
        k = ymax
        m=0.002
        #initialize vector
        p0 = (k, tau, m)
        popt, pcov = curve_fit(func, x[Indexs],y[Indexs], p0)
        
        
        if debug==1:
            
            # create a high resolution data set for the fitted waveform
            x2 = np.linspace(x[min(Indexs)], x[max(Indexs)], 7*10)
            y2 = func(x2, *popt)
            # plot the fitted curve and the real one
            plot(x2,y2)
            plot(x[Indexs],y[Indexs])
            
        return popt, pcov, tau    

# tests en cours    
if __name__ == "__main__":    
    #filename="/Users/Nicolas/Documents/Projets/TDO/post traitement/tracking freq/s120725d.sig"
    filename="/Users/nicolas/Desktop/s120725a.sig"
    print 'Importing file...'
    sig=loadfileGage(filename)
    print 'Done'
    # To reduce signals size
    sig2=sig[400000:4400000]
    t0=time.clock()
    #test differents algoritms
    print 'FFTGaussianfit...'
    Amp,Freq,FreqMax1,Sigma=FFTGaussianfit(sig)
    print 'Done in %f s' % (time.clock()-t0)
    #fints,intsb=ZeroCrossingTDO(sig)
    print "FFTZeropadTDO..."
    t0=time.clock()
    AmpMax,FreqMax=FFTZeropadTDO(sig2)
    print 'Done in %f s' % (time.clock()-t0)
    
    plt.figure()
    plt.plot(Freq)
    #plot(Freq-Sigma)
    #plot(Freq+Sigma)
    plt.plot(FreqMax1)
    #plt.plot(FreqMax)
    plt.show()
      